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ABSTRACT 

A previously-described code for constructing moving reversing layers (MRL) is improved by replacing a two-parameter 
model for g l (v), the radiative acceleration due to lines, with a flexible non-parametric description, thus allowing a 
greater degree of dynamical consistency to be achieved in modelling turbulent transonic flow in the outer atmospheric 
layers of O stars. With this new code, mass fluxes J are computed at fifty-seven points in (Tea, <?)-space. Specifically, 
J's are computed for all Lanz- Hubeny (2003) NLTE atmospheres with T cS (kK) £ (27.5,55) and log g(cm s" 2 ) < 4.5. 
Differences with widely-used mass-loss formulae are emphasized, and opportunities for differential spectroscopic tests 
identified. 
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1. Introduction 

In an earlier paper (Lucy 2007; L07), the moving reversing 
layer (MRL) theory of Lucy & Solomon (1970; LS70) was 
updated by incorporating an extensive line list and improv- 
ing the treatment of line formation. In addition, the models 
were, in effect, grafted onto the TLUSTY static NLTE O- 
star atmospheres of Lanz & Hubeny (2003) by imposing 
the TLUSTY emergent continuum flux distribution as the 
lower radiative boundary condition and by matching ion- 
ization fractions at T re 0.75T c ff. 

The motivation for reviving MRL theory was the con- 
flict between observed and predicted mass loss rates ($), 
which had led several spectroscopic groups to question the 
theory of radiatively-driven winds (e.g., Bouret et al. 2005; 
Fullerton et al. 2006). But, as emphasized in L07, the $'s 
being tested were not obtained by solving the equations 
governing the dynamics of radiatively-driven winds but are 
the values <&y derived by Vink et al. (2000) with a refined 
version of the semi-empirical Monte Carlo (MC) method 
of Abbott & Lucy (1985). Logically, therefore, the con- 
flict could arise from the specific assumptions of Vink et 
al. rather than from a failure of the radiative-driving mech- 
anism. 

In L07, MRL models were used to explore the sensitiv- 
ity of the eigenvalue J — Q/4irR 2 to Vt, the microturbulent 
velocity. Crucially, this parameter affects the flux irradiat- 
ing lines as matter is driven through the sonic point; and 
an increase in u t from 10 to 15 km/s was found to decrease 
J by re 0.3 dex. Thus, a physical effect was identified that, 
by reducing the predicted $'s, might partially contribute 
to resolving the conflict. 

More recently (Lucy 2010; L10), MRL theory was used 
to investigate individual stars. This was prompted by the 
work of Marcolino et al. (2009) on the weak wind problem, 
the major discrepancy for late-type O dwarfs between their 
$y's and observational estimates. To investigate this, a grid 



of 29 models was computed from which J's for particular 
stars could be obtained by interpolation. The results were 
encouraging: although the extremely low and uncertain J's 
estimated by Marcolino et al. were not matched, the pre- 
dicted J's were re 1.4 dex lower than the Jy's. Moreover, 
when J's were interpolated for the two strong-wind 04 stars 
analysed by Bouret et al.(2005), the results were consistent 
with the low values found when these authors took wind 
clumping into account. 

In the above investigation, the L07 code was deliber- 
ately not changed, thereby avoiding any suspicion that ad- 
justments were motivated by the observational data requir- 
ing explanation. But with the technique's usefulness thus 
demonstrated, a possibly significant flaw is now addressed, 
namely the local departures from dynamical consistency 
that result from the simple two-parameter representation of 
line driving. This further development is especially appro- 
priate since dynamical consistency in modelling transonic 
flow was identified in L10 as the key to accurate predic- 
tions of J and $. Accordingly, the primary purposes of this 
paper are first to describe how such improved models can 
be constructed and then to compute J's for all TLUSTY 
atmospheres relevant for Galactic main-sequence O stars. 

Throughout this paper J's are in units gm / s / cm 2 and 
$'s in M@/yr. 



2. Improved solution technique 

As in previous papers, transonic flow is assumed to be sta- 
tionary, isothermal, and plane-parallel. The equation of mo- 
tion can then be written as 



(v 2 



1 dv 
v dx 



"Seff 



(1) 
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Here a is the isothermal speed of sound, and g c s = g — g e — 
g l is the effective gravity, where g e = T e g and g l are the 
radiative accelerations due to electron- and line scatterings, 
respectively. 
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For given stellar parameters, we wish to find the solution 
of Eq. (1) such that the flow accelerates smoothly from 
sub- to supersonic velocities. This is achieved by finding 
the particular mass flux J that gives g c g = at the sonic 
point v = a, thus avoiding a singularity when integrating 
Eq. (1). 

When approximated by a MC estimator, g (v) is not 
analytic, and so solving Eq. (1) and determining its eigen- 
value J is not a conventional excercisc in integrating an 
ODE. Accordingly, in L07, a two-parameter formula was 
adopted for g l (v) that automatically gives g c e — at v = a, 
thus allowing the singularity-free stratification of the MRL 
to be obtained with a conventional integration of Eq. (1). 
The MC transfer calculation was then carried out in this 
stratified medium, resulting in estimates g l for each layer 
of the MRL. The challenge then was to find the values of J 
and of the parameters S and s that brought g l into optimal 
agreement with g (v; J,S, s). 

Because of this parametric approach, the solutions ob- 
tained had noticeable residuals Ag l = g l — g l - see Fig. 3 
in L10 - implying some uncertainty in the predicted J's. 
Although evidence was presented that J is moderately in- 
sensitive to departures from detailed local dynamical con- 
sistency, it is nevertheless desirable to eliminate this weak 
point in MRL theory. 

2.1. Non-parametric g l (v) 

To allow the g l, s to be accurately modelled by g l (v), this 
function is constrained to pass through the discrete set of 
points (gl, v^, with v\ < vi < ... < vi. The required contin- 
uous function g l (v) is then constructed as follows: by linear 
logarithmic interpolation between neighbouring points for 
v e (v\,vi); by setting g l (v) — g[ for v < V\; and by ex- 
trapolating the power law from the interval (7 — 1,/) for 
v > vj. 

The discrete representation extends from a small sub- 
sonic velocity v\ to a supersonic velocity vi > 2a, with 
spacing chosen to model the often sharply changing veloc- 
ity gradient as the sonic point is approached - see Fig.l 
in L07. The fc-th point is located at the sonic point - i.e., 
Vk = a - and the corresponding value of g\. is constrained 
to be <?* = g — g e , so that the regularity condition is again 
automatically satisfied. 

2.2. Stratification 

With the g l (v) thus defined, the MRL's stratification is 
obtained as described in Sect. 2.3 of L07, namely by two 
initial-value integrations of Eq. (1) starting at v = a, one 
for v < a and one for v > a. 

To avoid a singularity at v = a, the initial velocity gra- 
dient must be such that 



dv 
i — 

dx 
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Now, if g (v) were an analytic function, the logarithmic 
derivative in Eq. (2) would be the same for both inward 
and outward integrations. But the adopted piecewise-lincar 
segmented representation of logg 1 is not analytic: although 
continuous, its derivative is in general discontinuous at 
v = v i and thus may be so at v k = a. Accordingly, for the 



inward and outward integrations, the logarithmic deriva- 
tives are the slopes ofthe ( k — 1 , k ) and (k, k + 1 ) segments , 
respectively. 

The discontinuities in the derivatives of g (v) allow the 
representation to approximate curvature in (log g — logv)- 
plots. Of course, as / — > oo, unlimited accuracy can be 
achieved, and the discontinuities then — »■ 0. 

2.3. An example 

Model 44005375, with parameters T off = A0,000K,log g = 
3.75 and v t — 10km/ s, illustrates the improved technique. 

Fig.l shows the first steps in the search for J. The start- 
ing values for g\^ k are obtained from the two-parameter 
formula - see Eq. (1) in L10 - with S — 0.5 and s = 1.5 -i.e., 
a broken power law, with the switch to a rapidly increasing 
g (v) occurring at Mach number m — v/a = 0.63. With 
g l (v) thus fixed, several models were computed with vary- 
ing J in order to locate the root of Qi^(J) = - see Sect. 2. 3 
and Fig. 2 in L10. The result, J = —5.64 dex, is then such 
that, as matter accelerates from mi = 0.5 to m 2 = 2.0, 
the work done by the gradients of gas and radiation pres- 
sures accounts for the gain in mechanical energy per gm. 
Nevertheless, the non-vanishing residuals Ag l imply that 
the I — 1 values g\. k require adjustment. 

Ideally, corrections to gj^ k should be derived from 

the Ag l with an algorithm analogous to the temperature- 
correction procedures in stellar atmosphere theory. But 
here a trial-and-error procedure is followed based on in- 
spection of plots such as Fig.l. Thus, Fig.l shows that a 
steeper slope than s = 1.5 is required for m > 1 and that 
the g\ should be increased by w 0.1 dex for m e (0.1,0.4). 

With the g\, k thus adjusted, a new sequence of models 
is computed, the modified root of = derived, and an 
updated version of Fig.l plotted. This iterative procedure 
is continued until a satisfactory degree of convergence is 
achieved. In this case, the final model has J = —5.72 dex 
and is plotted in Fig. 2. The iteratively-corrected function 
g l (v) now agrees closely with the g l, s. A dynamically con- 
sistent model of turbulent transonic flow has therefore been 
constructed. 



3. Microturbulence 

As demonstrated in L07, the J's predicted by MRL theory 
are sensitive to v t , which must therefore be included with 
T e g and g when comparing with observational data. 

Given the importance of this sensitivity, the improved 
code is now applied to the model £400g375 in order to check 
and extend the analysis of L07. As in that investigation, 
when Vt is varied, the b- values and incident flux distribution 
are from the TLUSTY model with v t = lOkm/s. 

Solutions ranging from v t — - i.e., pure thermal broad- 
ening in the lines' Doppler cores - to v t — 20km/ s, corre- 
sponding to near sonic turbulence, are plotted in Fig. 3, to- 
gether with the two solutions from L07 at 6.7 and 10km/ s. 
The sensitivity to v t is confirmed. 

From the new results between 6.7 and 15km/ s, the log- 
arithmic slope at 10km /s is 



d log J 
d logv t 



= -1.46 



(3) 
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Fig. 1. Search for dynamical consistency. The open circles are 
the initial discrete representation (g\,Vi) for model i400(/375, 
and the connecting solid line is the resulting continuous func- 
tion g l (v). The MC estimates g l are plotted as filled circles and 
correspond to the mass flux J = —5.64 dex that gives Qi,2 = 0. 
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Fig. 2. Search for dynamical consistency. Symbols and model 
parameters as in Fig.l. The interatively-corrected representation 
(gl,Vi) is shown together with the resulting MC estimates g l 
when J = —5.72 dex. 



This value is recommended for propagation-of-error calcula- 
tions oi<Ji gj for stars with u t 's comparable to the canonical 
v t = 10km/ s. 

Fig. 3 shows that, in principle, the parameter vt could 
bring about differences in J by ^ 1 dex at fixed T e s,g. 
But as this seems not to happen for real stars, the mecha- 
nism exciting and maintaining microturbulence is presum- 
ably preventing the independent variation of this parame- 
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Fig. 3. Sensitivity of the eigenvalue J to microturbulent velocity 
v t for model £400*7375. Filled circles are solutions obtained as 
described in Sect. 2; open circles are solutions from L07. The 
thermal speed of Fe ions is indicated. 



ter, resulting perhaps in a functional dependence of v t on 
Tes,g. 

In view of the strong damping expected for turbulence 
with v t /a ~ 0.5, work done by radiative forces is almost cer- 
tainly required for its maintenance, in which case the phe- 
nomenon of microturbulence is not unrelated to radiatively- 
driven outflows. In fact, Fig. 3 suggests a direct, causal re- 
lationship as follows: Given the observational evidence that 
wind-clumping occurs shortly after the sonic point (Bouret 
et al. 2005), some clumps may well lose their net outward 
driving and thus fall back into the photosphere (e.g., Howk 
et al. 2000) where their dissipated kinetic energy could ex- 
cite and maintain local turbulence. If the fall-back fraction 
decreases with decreasing J, a feedback loop operates, so 
that for a given star only one pairing (J,Vt) is possible. 

This conjectured coupling of J and v t is supported by 
the difficulty of obtaining satisfactory high mass flux MRL's 
when v t — > 0. The J's for v t ^ 5fcm/s in Fig. 3 are esti- 
mated by achieving dynamical consistency only for v/a < 1 
since the extension to v/a ^ 2 is not possible because g e g 
becomes positive - see Fig. 4. Admittedly, this breakdown 
of the iterative procedure is code-specific, resulting from 
choosing v rather than height x as the independent vari- 
able. Nevertheless, even if the code were reconfigured to 
allow dynamical consistency to be extended to supersonic 
velocities, the resulting non-monotonic velocity law implies 
a density inversion at v/a ~ 2 — 3. Such a stratification is 
surely unstable and thus a possible origin of clumping and 
of infalling blobs. 

3.1. A spectroscopic test 

This sensitivity of J to vt is in stark contrast to CAK the- 
ory (Castor, Abbottfe Klein 1975), which posits that the 
properties of radiatively-driven winds can be derived on 
the basis of the Sobolev approximation. In this approxima- 
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Fig. 4. Failed search for dynamical consistency when vt = for 
model t400p375. Symbols as in Fig.l. 

tion, g l (v) is independent of the lines' absorption profiles 
and therefore independent of v t 's contribution to the width 
of the Doppler core. 

In principle, this can be tested observationally, and is 
best done differentially. Thus stars differing in vt but not 
widely separated in (T cff , g)-space should be subjected to 
identical observing and diagnostic procedures to see if A$ 
can be understood without (CAK), or only with (MRL), a 
contribition from Av t . 

4. Computed mass fluxes 

In this section, the improved technique of Sect. 2 is used 
to recalculate the 29 models in Table 1 of L10 and then 
to add a further 28 models in order to provide a rather 
complete coverage of (T c ff , g)-space for H-burning O stars. 
As before, the models' composition is solar with Nhc/Nh = 
O.f (Grevess & Sauval 1998), the included metal ions are as 
in Table 1 of Lanz & Hubeny (2003), and v t = 10km/ s. 

4.1. O-star grid 

Mass fluxes J for 57 models are given in Table 1. The grid 
is determined by the availability of TLUSTY atmospheres 
(Lanz & Hubeny 2003) and is complete for their models 
with logg < 4.5. Because O stars on the ZAMS have logg « 
4.2 - see Fig. 4 in L10, J's have not been computed for 
log g = 4.75. 

The coverage provided by Table 1 allows J's to be deter- 
mined for all Galactic O stars by interpolation - or a slight 
extrapolation in a few cases. If, in addition to T c g, g and v tl 
a star's distance is known, its radius can be computed and 
therefore also $ = AnR J. Given the detailed diagnostic 
modelling of numerous O stars in recent years, a critical 
evaluation of MRL theory may be possible with existing 
data. This is not attempted here and is, in any case, best 
carried out by investigators familiar with the uncertainties 
of analysing circumstellar spectra. 
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J — J* defines the boundary of the weak-wind domain - see 
L10. 
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Fig. 6. Mass flux J as a function of g for T e g = 27.5(2. 5)55. OkK. 
The data are from Table 1. 



The data in Table 1 is also relevant for investigations of 
stellar evolution with mass loss, for computing the latitude 
dependence of mass loss for rapidly- rotating stars, and for 
calculating its radial dependence for accretion disks. 

The dependence of J(T e s,g;vt) on T c ff and g is shown 
in Figs. 5 and 6. The most striking feature is the departure 
from the expected monotonic increase of J with increasing 
T e s that occurs when T e g ^ 30,000-fT and logg <^ 3.9. As 
discussed in L10, this prediction of MRL theory offers at 
least a partial explanation of the weak-wind phenomenon. 
Interestingly, the minimum at T e s « 30, 000K is deeper for 
g > gzAMS 
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Table 1. Computed mass fluxes J. 
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4.2. Accuracy 

The J's in Table 1 have several sources of uncertainty. The 
first originates from MC sampling errors - see Appendix 
A. The g l, a sampling errors propagate via Qi t 2 into errors 
in J when this quantity is determined by locating the in- 
tercept Qi^(J) — - see Fig.2 in L10. Thus, for the fi- 
nal M00c/375 model in Fig.2, the least squares solution is 
log J = —5.715 ± 0.014. This estimate of ai og j shows that 
MC noise is inconsequential since other errors are surely far 
greater. 

The uncertainty ai og j — > as the number of MC quanta 
N — > oo, but J would still be subject to error because of 
residual departures from exact dynamical consistency if the 
continuous function g l (v) derives from only a finite number 
of points {g\,Vi). 

This second source of uncertainty can be estimated as 
follows: Given that the non-parametric representation is a 
marked improvement over the previous model, the average 
\AlogJ\ between the models treated both here and in L10 
is a good estimate of the typical error of the previous J's 
and, at the same time, a conservative error estimate for the 
J's in Table 1. The 29 J's differ on average by only 0.10 
dex, with the largest difference being 0.33 dex. 

A third source of uncertainty is the cumulative effect of 
errors in abundances, input physics and line- formation the- 
ory, many of which propagate from the TLUSTY models. 
A plausible guess is that these errors should rarely exceed 
0.2 dex. 

A fourth source of uncertainty is the derivation of J 
from a plane-parallel treatment of transonic flow, with 



back-scattering from v > 5a neglected. This is investigated 
in Appendix B and found to well-justified. 

A fifth and probably dominant source of uncertainty 
is the sensitivity to the throttling effect of turbulent line 
broadening (L07 and Sect. 3 above). Given our ignorance as 
to the source and nature of photospheric turbulence, this 
has perforce been investigated in the microturbulent limit, 
and strong sensitivity is found. Thus, from Fig. 3, we see 
that, with v t — lOkm/s, photospheric turbulence reduces 
J by ~ 1.3 dex from its value for laminar flow - i.e., pure 
thermal broadening. Within the context of the microturbu- 
lent model, an error of ±2fcm/s at v t ~ 10km/ sec implies 
o~iogj ~ 0.13 according to Eq. (3). 

Given these uncertainties, tests of MRL theory, either 
spectroscopically or via stellar evolution calculations, might 
reasonably allow for an uncertainty of ±0.2 dex in the J's 
given in Table 1. But if a test reveals systematic discrepan- 
cies > 0.4 dex, a contradiction can be claimed. 

5. Comparisons with mass-loss formulae 

In this section, the MRL mass fluxes are compared to 
widely-used mass-loss formulae. The aim here is not a com- 
prehensive discussion but to call attention to opportunities 
for differential testing, as already discussed in Sect. 3.1. 

5.1. CAK 

In their recent discussion of WNH stars, Smith & Conti 
(2008) sketch an evolutionary scenario based on an O-star 
mass-loss formula extracted from CAK theory. For fixed 
T e g, their formula gives J's dependence on g = g e /T e as 

where Jo is the mass flux when T e = 0.5. Eq. (4) predicts 
that J — > oo as T e — > 1 - i.e., as the Eddington limit is 
approached - and this is basic to their claim that a feed- 
back process results in runaway mass loss late in the core-H 
burning phase of massive O stars. 

In Fig. 7, the behaviour of Jl as g — > g e is compared 
to Eq.(4) when J is chosen to match Jl at log g = 4. 
This comparison is carried out for T e g = 30, 35 and 
40kK, for which the lowest g TLUSTY models have T e = 
0.52, 0.54 and 0.52, respectively. Remarkably, MRL theory 
predicts greatly enhanced mass loss when a star is still well- 
detatched from its Eddington limit. In each of the plotted 
sequences, the lowest g model is offset from the scaled CAK 
formula by A J ~ 1 dex. 

If the Smith-Conti scenario were supported by actual 
stellar evolution tracks and accurately fitted several well- 
observed WNH stars, then MRL theory could be immedi- 
ately dismissed as over-predicting J for evolved O stars. But 
neither of these circumstances obtains, and this enhanced 
mass loss already at T e m 0.5 will likely prove favourable 
to their scenario. Clearly, further work is required on the 
evolution of mass-losing massive stars. 

In addition to its relevance for evolutionary scenarios, 
Fig. 7 suggests another differential test to distinguish MRL 
and CAK theories. In this case, pairs of stars with closely 
similar Tiff's but markedly different g's should be observed 
and analysed identically to see which ^-dependence in Fig. 
7 is favoured. 
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Fig. 7. Comparison of MRL mass fluxes J with scaled CAK 
formulae for T e g = 30, 35 and 40kK. The dashed vertical lines 
are the MRL-CAK offsets discussed in the text. 
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Fig. 8. Mass-loss rates for O-stars on the zero-age main sequence 
(ZAMS) as a function of T e g, with masses in solar units indi- 
cated. Predictions of MRL theory (<&l) are plotted as well as 
values (3>v) derived with the formula of Vink et al. (2000). 



5.2. Vink et al. 

A model from Table 1 cannot be directly compared to the 
Vink et al. formula since the latter requires three fundamen- 
tal stellar parameters and Voc . Accordingly, the comparison 
is carried out for the ZAMS models of Pols et al. (1998) 
with Z = 0.02. Each point on the ZAMS gives M,R and L 
which, with the additional assumption that Voo/v esc = 2.6 
(Lamers et al. 1995), allows $y to be computed from Eq. 
(12) of Vink et al. (2000). This can then be compared to 
$l = 47ri? 2 x J(T e g,g). Here J is obtained by simple bi- 
variate interpolation (Abramowitz & Stegun 1965) from the 
four surrounding entries in Table 1, taking the independent 
variables to be log T e g and log g. 

The two predictions for the ZAMS are plotted in Fig. 
8. Throughout the entire range A$ = <Ev - $ L > 0. At 
T c g — 27, 650K, the offset is 0.40 dex, and this increases 
steeply to a maximum of 1.24 dex at 30, 150i^T. Thereafter, 
A$ decreases - not quite monotonically - to reach a barely 
significant 0.27 dex at T cff = 50, 300^. 

The huge difference at T off w 30, OOOif allows MRL the- 
ory to partly explain the weak wind phenomenon, which 
arose when the Vink et al. predictions were compared to 
mass loss estimates for late-type O dwarfs. In fact, the rel- 
evant diagnostic analyses (Marcolino et al. 2009) appear 
still to require $l to be reduced by « 0.8 dex - see L10. 

According to Fig. 8, MRL theory predicts that a massive 
star's initial $ is markedly less than $y. However, when its 
expanding radius has reduced g by w 1 dex, the situation 
reverses- see Fig. 7. Thus, for example, if a mass-losing star 
reaches the point (40, OOO-ft', 3.5) in (T c e, log <?)-space with 
M/M Q = 80, then $y = -4.75 dex, but $ L = -4.59 dex, 
a factor 1.4 larger. To put this in context, the 120.Mq track 
of Limongi & Chieffi (2006) computed with Vink et al. mass 
loss has logg = 3.56 and 3.31 with mass 93.6 and 56.6.Mq, 
respectively, at its two crossings of T c g = 40,000if during 
core H-burning. 



From the trends evident in Fig. 7, the MRL mass-loss 
enhancement will be even greater for g's smaller than given 
in Table 1. Moreover, the Limongi-Chieffi track suggests 
that this part of parameter space may well be accessed by 
real stars. Accordingly, Table 1 needs to be extended to 
lower g's. 

6. Conclusion 

Motivated by the (partial) successes of MRL theory in re- 
producing the reduced $'s found by spectroscopists for O 
stars when wind-clumping is taken into account, the aim of 
this paper has been to complete the coverage of the relevant 
(Teg , (;)-domain provided by the TLUSTY atmospheres. To 
this end, the MRL code has been improved by adopting a 
non-parametric description of g l (v) 7 which has allowed a 
greater degree of dynamical consistency to be achieved for 
turbulent transonic flow, and thus more accurate J's. The 
results of this effort are the 57 values of log J in Table 1. 

Interpolation in Table 1 allows J to be derived for any O 
star with measured T e g and g. Moreover, a propagation-of- 
error calculation gives ai og j if the standard errors of these 
two parameters and of Vt have been determined. This can 
be done for stars individually and independently. However, 
given the difficulties of diagnostic analyses, with the result- 
ing possibility of systematic errors, there is merit in per- 
forming differential tests as suggested in Sects. 3.1 and 5.1. 
If v t is a slowly-varying function of T e g and g, the first test 
may be indecisive. But the proposed test between the MRL 
and CAK theories suggested by Fig. 7 is feasible. Moreover, 
this test is fundamental for stellar wind theory since it di- 
rectly concerns the question: where in the outflow is the 
mass-loss rate determined? 

With regard to this question, the partial success of MRL 
theory in explaining the weak wind phenomenon supports 
the LS70 argument that J and therefore <f> is determined by 
the regularity condition at the sonic point. This would be 
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Fig. A.l. Convergence and accuracy test of MC estimators. For 
model t500g400, the mean absolute fractional error e of g l for 
0.5 < v/a < 2 is plotted against N, the number of energy pack- 
ets. The filled and open circles refer to estimators A and B, 
repectively. 



decisive if the observed $'s of the Marcolino et al. (2009) 
stars were convincingly revised into agreement with the pre- 
dictions of MRL theory. But since only weak C IV absorp- 
tion is observed, there is little or no observational basis for 
such improved estimates. Other tests should therefore be 
carried out. 
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Appendix A: Precision of the estimator g l 

The g l are derived from a MC simulation using estimator A, the 
summation over pathlengths given in Eq. (10) of L07. This is expected 
to be more accurate than estimator B, the summation of momentum 
tansfers from energy packets to matter at the discrete line-scattering 
events. 

A test of the accuracy and convergence of A and B has been 
carried out for model <500<?400. In this test, the exact g l is taken to 
be g l given by A when N, the number of MC quanta, is 512 X 10 6 . 
Given this 'exact' g l , the fractional errors of g l can be computed at 
smaller N for A and at all N for B. 

The mean absolute fractional errors for 0.5 < v/a < 2 arc plot- 
ted against N in Fig. A.l. As expected, errors for both A and B arc 
<x l/vN. Also as expected, A is the more accurate. From the plotted 
least squares fits, e\ = 0.488 X eg. Thus, to achieve the same accu- 
racy, B would have required N to be increased by the factor 4.2. The 
saving of computer time with A was essential in carrying out the huge 
modelling effort required to produce Table 1. 

A typical simulation has N = 40 X 10 6 and so, according to Fig. 
A.l., the typical error of g l with A is 0.016 dex. But this is specific 
to model 45003400. At cooler T e fj's, an increasing fraction of packets 
propagate through the MRL without undergoing line scatterings, so 
that, for fixed N, the sampling error of g l increases. Thus, for model 
t325<?400, the above error increases to 0.20 dex. 



Appendix B: Reflection probability 

When a MC quantum crosses the upper boundary of the MRL at 
height x\ where v/a = 5, it permanently exits the computational 
domain. Thus, its dimcnsionlcss reflection probability pi = 0, an as- 
sumption justified if the true pi < 1. 

An estimate of pi can be derived following the method of Abbott 
& Lucy (1985; Sect III). We assume that radiative driving dominates 
other mechanisms beyond x\ and accelerates the wind to terminal ve- 
locity Voa . On this assumption, increases in the mechanical luminosity 



C(r) = $ ( 



1 



GM 



(B.l) 



are accounted for by a matching decrease in the radiative luminosity 
L(r). Thus, in the spherical shell (r, r + dr), the 0(v/c) difference 
between the rates at which matter absorbs dA and emits dS radiant 
energy is 

dA-d£ = dC = $ (v— + g) dr (B.2) 
dr 

If we now assume that the electron- and linc-scattcrings responsi- 
ble for the energy transfer absorb from a radially-streaming radiation 
field and emit isotropically, then 



d€ 



dC 



(B.3) 



Since this energy is radiated isotropically, the fraction propagating 
back through r = n - i.e., back into the MRL - is 



V 1 - (*/*l) 2 



(B.4) 



provided that no interactions intervene. Here z = R/r, where R is the 
photosphcric radius. 

Combining the above, we find that the luminosity of inwardly- 
propagating radiation at r\ is 



J n 



c dC , 
wj — dr 



v dr 



(B.5) 



The reflection probability is then pi 
the luminosity the wind-free star. 

In evaluating pi, we assume L~[ 
wind obeys a /3-vclocity law with fS = 



L 1 /(L* + L x ), where L* is 



-C L* and that the supersonic 
1. The result is 



Pi 



t Jo 



W! [ 1 + ■ 



dz 



(B.6) 



Here 4>| = L^/cVoc (Cassinelli & Castor 1973) and r\ = (u e sc/foo) 2 , 
where v eac is the escape velocity from r = R. 

Values of pi have been computed for the ZAMS models in Fig. 8. 
Thus, at M = 30Mq with v^/v^ = 2.6, the ratio $ L /$i = 0.044. 
Then, for v\ = 5a, z\ = 0.969, and we find pi = 0.0071. The neglect of 
back-scattering from the exterior wind is therefore justified. Moreover, 
since 1 — z\ -C 1, the assumption of plane-parallel geometry is also 
justified. 

But for very massive stars, the assumption pi = is less valid. 
Thus, for M = 8OA40, *z,/<I>t = 0.384, z\ = 0.972, and we find that 
pi = 0.064. Accordingly, when theory and observation agree to ^ 0.1 
dex, further progress will require an improved treatment of transonic 
flow. 
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